#include <admodel.h>

  extern "C"  {
    void ad_boundf(int i);
  }
#include <tm_M.htp>

model_data::model_data(int argc,char * argv[]) : ad_comm(argc,argv)
{
  styr.allocate("styr");
  endyr.allocate("endyr");
 n_model_years = endyr-styr+1;
  n_intrvls_per_yr.allocate("n_intrvls_per_yr");
  fishery_sample_prob.allocate(1,n_model_years,1,n_intrvls_per_yr,"fishery_sample_prob");
  survey_sample_prob.allocate(1,n_model_years,1,n_intrvls_per_yr,"survey_sample_prob");
  survey_effort.allocate(1,n_model_years,1,n_intrvls_per_yr,"survey_effort");
  fishery_sample_prob_vec.allocate(1,n_model_years*n_intrvls_per_yr);
  survey_sample_prob_vec.allocate(1,n_model_years*n_intrvls_per_yr);
  survey_effort_vec.allocate(1,n_model_years*n_intrvls_per_yr);
  max_age.allocate("max_age");
 first_cohort = styr-max_age+1;  
 last_cohort = endyr;  
  fishery_n_obs.allocate("fishery_n_obs");
  fishery_dat.allocate(1,fishery_n_obs,1,4,"fishery_dat");
  survey_n_obs.allocate("survey_n_obs");
  survey_dat.allocate(1,survey_n_obs,1,4,"survey_dat");
  fishdat_tend.allocate(1,fishery_n_obs);
  fishdat_age_mid.allocate(1,fishery_n_obs);
  fishdat_catch.allocate(1,fishery_n_obs);
  fishdat_cohort.allocate(1,fishery_n_obs);
  fishdat_cohort_index.allocate(1,fishery_n_obs);
  fishdat_cohort_intrvl_index.allocate(1,fishery_n_obs);
  fishdat_year.allocate(1,fishery_n_obs);
  fishdat_month.allocate(1,fishery_n_obs);
  survey_tend.allocate(1,survey_n_obs);
  survey_age_mid.allocate(1,survey_n_obs);
  survey_catch.allocate(1,survey_n_obs);
  survey_cohort.allocate(1,survey_n_obs);
  survey_cohort_index.allocate(1,survey_n_obs);
  survey_cohort_intrvl_index.allocate(1,fishery_n_obs);
  survey_year.allocate(1,survey_n_obs);
  survey_month.allocate(1,survey_n_obs);
 n_cohorts = n_model_years+max_age-1; 
 cout << n_cohorts << " " << first_cohort << " " << last_cohort << " " << last_cohort - first_cohort + 1 << endl;
  n_intrvls_cohort.allocate(1,n_cohorts);
  beg_intrvl_cohort.allocate(1,n_cohorts);
  end_intrvl_cohort.allocate(1,n_cohorts);
  n_model_intrvls = (endyr-styr+1)*n_intrvls_per_yr;
  // Get data into vectors (from matrix
  fishdat_tend    = column(fishery_dat,1);
  fishdat_age_mid = column(fishery_dat,2);
  fishdat_catch     = column(fishery_dat,3);
  fishdat_cohort     = column(fishery_dat,4);
  for (int i=1;i<=fishery_n_obs;i++)
  {
    fishdat_year(i) = int(fishdat_tend(i) - .001);    
    fishdat_month(i) = n_intrvls_per_yr*(fishdat_tend(i) - fishdat_year(i)+.001);    
  }
  survey_sample_prob_vec.initialize();
  fishery_sample_prob_vec.initialize();
  survey_effort_vec.initialize();
  for (int i=1;i<=n_model_years;i++)
  {
  /*
   	cout << n_model_years << endl;
  	cout << i <<endl;
  	cout << (i-1)*n_intrvls_per_yr+1 << " " << i*n_intrvls_per_yr << endl;
  	cout << survey_sample_prob(i) << endl;
  	cout << survey_sample_prob_vec((i-1)*n_intrvls_per_yr+1,i*n_intrvls_per_yr) << endl;
  */
  	for(int j = 1; j <=n_intrvls_per_yr;j++)
  	{
      survey_sample_prob_vec((i-1)*n_intrvls_per_yr+j) = survey_sample_prob(i,j);
  	  //cout << survey_sample_prob_vec((i-1)*n_intrvls_per_yr+j) << endl;
      survey_effort_vec((i-1)*n_intrvls_per_yr+j) = survey_effort(i,j);
      fishery_sample_prob_vec((i-1)*n_intrvls_per_yr+j) = fishery_sample_prob(i,j);
    }
  }
  survey_tend    = column(survey_dat,1);
  survey_age_mid = column(survey_dat,2);
  survey_catch     = column(survey_dat,3);
  survey_cohort     = column(survey_dat,4);
  for (int i=1;i<=survey_n_obs;i++)
  {
    survey_year(i) = int(survey_tend(i) - .001);    
    survey_month(i) = n_intrvls_per_yr*(survey_tend(i) - survey_year(i)+.001);    
  }
  // cout <<fishdat_year<<endl;
  //cout <<survey_year<<endl;
  //cout <<fishdat_month<<endl;
  //cout <<survey_month<<endl;
  // Vector of interval definitions
  
  // Compute the number of observations by 
  //nages = max_age;
  //maxK  = 1+endyr-styr+max_age;
  //minK  = 1;
  Min=0.2;
  //offset = 0.0;
  end_intrvl_cohort.initialize();
  beg_intrvl_cohort.initialize();
  n_intrvls_cohort.initialize();
  for(int k=1;k<=max_age;k++)
  {
    beg_intrvl_cohort(k) = 1;
    end_intrvl_cohort(k) = k*n_intrvls_per_yr;
    n_intrvls_cohort(k)  = end_intrvl_cohort(k) - beg_intrvl_cohort(k) + 1;
  }
  for(int k=max_age;k<=n_cohorts;k++)
  {
    beg_intrvl_cohort(k) = (k-max_age)*n_intrvls_per_yr + 1 ;
    if(max_age*n_intrvls_per_yr < n_model_intrvls - beg_intrvl_cohort(k)+1.)
      end_intrvl_cohort(k) = beg_intrvl_cohort(k) + max_age*n_intrvls_per_yr;
    else
      end_intrvl_cohort(k) = n_model_intrvls;
    // Number of intervals for each cohort
    n_intrvls_cohort(k)  = end_intrvl_cohort(k) - beg_intrvl_cohort(k) + 1;
  } 
  // cout<<n_intrvls_cohort<<endl;
  // cout<<beg_intrvl_cohort<<endl;
  // cout<<end_intrvl_cohort<<endl;
  // cout<<n_cohorts.indexmax()<<endl;
  Ages.allocate(1,max_age*n_intrvls_per_yr);
  ages_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort);
  cal_yr_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort);
  cal_intrvl_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort);
  fishery_catch_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort);
  survey_catch_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort);
  total_fishery_catch_cohort.allocate(1,n_cohorts);
  total_survey_catch_cohort.allocate(1,n_cohorts);
  for (int i=1;i<=n_intrvls_per_yr*max_age;i++)
    Ages(i) = double((i-.5)/double(n_intrvls_per_yr));
  // cout << Ages<<endl;exit(1);
  fishdat_year.initialize();
  survey_year.initialize();
  total_survey_catch_cohort.initialize();
  survey_catch_cohort.initialize();
  total_fishery_catch_cohort.initialize();
  fishery_catch_cohort.initialize();
  int ii;
  for (int i=1;i<=survey_n_obs;i++)
  {
    //int k  = survey_cohort(i)-styr+max_age;
   //cout <<"i: " << i << " of " << survey_n_obs <<endl;
   survey_cohort_index(i) = survey_cohort(i)-styr+max_age;
   if (survey_cohort_index(i) <= max_age)
      ii = n_intrvls_per_yr*(survey_tend(i) - styr +.001);
    else
      ii = n_intrvls_per_yr*(survey_tend(i) - survey_cohort(i)+.001);
    // Dindex_S(k,ii) = i;
    survey_cohort_intrvl_index(i) = ii;
    survey_catch_cohort(survey_cohort_index(i),ii) += survey_catch(i);
    total_survey_catch_cohort(survey_cohort_index(i))    += survey_catch(i);
  }
  for (int i=1;i<=fishery_n_obs;i++)
  {
    //int k  = fishdat_cohort(i)-styr+max_age;
    fishdat_cohort_index(i) = fishdat_cohort(i)-styr+max_age;
    if (fishdat_cohort_index(i) <= max_age)
      ii = n_intrvls_per_yr*(fishdat_tend(i) - styr +.001);
    else
      ii = n_intrvls_per_yr*(fishdat_tend(i) - fishdat_cohort(i)+.001);
    fishdat_cohort_intrvl_index(i) = ii;
    fishery_catch_cohort(fishdat_cohort_index(i),ii) += fishdat_catch(i);
    total_fishery_catch_cohort(fishdat_cohort_index(i)) += fishdat_catch(i);
  }
  for (int k=1;k<=n_cohorts;k++)
  {
   for (int i=1;i<=n_intrvls_cohort(k);i++)
    {
      if (k<=max_age)
      {
        // Indexes age for the cohort by interval
        ages_cohort(k,i) = n_intrvls_per_yr*(max_age - k) + i; 
        cal_yr_cohort(k,i) = styr + int((double(i)-.5)/double(n_intrvls_per_yr));
      }
      else
      {
        ages_cohort(k,i) =  i ;
        cal_yr_cohort(k,i) = k-max_age +styr + int((double(i)-.5)/double(n_intrvls_per_yr));
      }
      cal_intrvl_cohort(k,i) = n_intrvls_per_yr*(cal_yr_cohort(k,1) - styr) + i;
    }
  }
  // Print out the Catch by cohort
  // cout << cal_intrvl_cohort <<endl;exit(1);
  //cout <<"Hear"<<endl;
  //exit(1);
  N0_phase.allocate("N0_phase");
  F_phase.allocate("F_phase");
  q_phase.allocate("q_phase");
  fishery_sel_phase.allocate("fishery_sel_phase");
  survey_sel_phase.allocate("survey_sel_phase");
  M_phase.allocate("M_phase");
}

void model_parameters::initializationfunction(void)
{
  log_N0.set_initial_value(13.52393);
  log_F.set_initial_value(-1.203973);
  fishery_sel50.set_initial_value(5.);
  fishery_selslope.set_initial_value(1.);
  survey_sel50.set_initial_value(2.);
  survey_selslope.set_initial_value(5.);
  log_q.set_initial_value(-4.199705);
  log_M.set_initial_value(-1.609438);
}

model_parameters::model_parameters(int sz,int argc,char * argv[]) : 
 model_data(argc,argv) , function_minimizer(sz)
{
  initializationfunction();
  log_N0.allocate(1,n_cohorts,N0_phase,"log_N0");
  log_F.allocate(1,n_model_years,F_phase,"log_F");
  log_q.allocate(q_phase,"log_q");
  F_fishery_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort,"F_fishery_cohort");
  #ifndef NO_AD_INITIALIZE
    F_fishery_cohort.initialize();
  #endif
  F_survey_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort,"F_survey_cohort");
  #ifndef NO_AD_INITIALIZE
    F_survey_cohort.initialize();
  #endif
  Z_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort,"Z_cohort");
  #ifndef NO_AD_INITIALIZE
    Z_cohort.initialize();
  #endif
  survival_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort,"survival_cohort");
  #ifndef NO_AD_INITIALIZE
    survival_cohort.initialize();
  #endif
  cumulative_survival_cohort.allocate(1,n_cohorts,1,n_intrvls_cohort,"cumulative_survival_cohort");
  #ifndef NO_AD_INITIALIZE
    cumulative_survival_cohort.initialize();
  #endif
  cumulative_nat_mort_cohort.allocate(1,n_cohorts,"cumulative_nat_mort_cohort");
  #ifndef NO_AD_INITIALIZE
    cumulative_nat_mort_cohort.initialize();
  #endif
  F.allocate(1,n_model_years,"F");
  #ifndef NO_AD_INITIALIZE
    F.initialize();
  #endif
  q.allocate("q");
  #ifndef NO_AD_INITIALIZE
  q.initialize();
  #endif
  fishery_sel50.allocate(fishery_sel_phase,"fishery_sel50");
  fishery_selslope.allocate(fishery_sel_phase,"fishery_selslope");
  survey_sel50.allocate(survey_sel_phase,"survey_sel50");
  survey_selslope.allocate(survey_sel_phase,"survey_selslope");
  fishery_select.allocate(1,max_age*n_intrvls_per_yr,"fishery_select");
  #ifndef NO_AD_INITIALIZE
    fishery_select.initialize();
  #endif
  survey_select.allocate(1,max_age*n_intrvls_per_yr,"survey_select");
  #ifndef NO_AD_INITIALIZE
    survey_select.initialize();
  #endif
  N0.allocate(1,n_cohorts,"N0");
  #ifndef NO_AD_INITIALIZE
    N0.initialize();
  #endif
  log_M.allocate(M_phase,"log_M");
  M.allocate("M");
  #ifndef NO_AD_INITIALIZE
  M.initialize();
  #endif
  log_exploit.allocate("log_exploit");
  #ifndef NO_AD_INITIALIZE
  log_exploit.initialize();
  #endif
  ll_i.allocate("ll_i");
  #ifndef NO_AD_INITIALIZE
  ll_i.initialize();
  #endif
  ll.allocate("ll");
  #ifndef NO_AD_INITIALIZE
  ll.initialize();
  #endif
  nll.allocate("nll");
}

void model_parameters::userfunction(void)
{
  F = mfexp(log_F);
  q = mfexp(log_q);
  M = mfexp(log_M);
  F_fishery_cohort.initialize();
  // cumulative_survival_cohort = 1.;
  // Compute selectivity (over all monthly ages?)
  for (int i=1; i<=max_age*n_intrvls_per_yr; i++)
  {
    fishery_select(i) = 1./(1+mfexp(-fishery_selslope*(Ages(i) - fishery_sel50)));
    survey_select(i) = 1./(1+mfexp(-survey_selslope*(Ages(i) - survey_sel50)));
  }
  for (int k=1; k<=n_cohorts; k++)
  {
    //int iy   = cal_yr_cohort(k,1); // assign temporary if using multiple times for overhead reasons
    int iy = cal_yr_cohort(k,1)- styr + 1;
    int i_intrvl = cal_intrvl_cohort(k,1);
    int iage = ages_cohort(k,1);
    F_fishery_cohort(k,1) = fishery_select(iage) * F(iy) * fishery_sample_prob_vec(i_intrvl);
    F_survey_cohort(k,1) = survey_select(iage) * q * survey_sample_prob_vec(i_intrvl) * survey_effort_vec(i_intrvl);
    Z_cohort(k,1) = F_fishery_cohort(k,1) + F_survey_cohort(k,1) + M;
    survival_cohort(k,1) = mfexp(-Z_cohort(k,1)/double(n_intrvls_per_yr));
    cumulative_survival_cohort(k,1) = survival_cohort(k,1);
    cumulative_nat_mort_cohort(k) = (1. - survival_cohort(k,1))*M/Z_cohort(k,1);
    for (int i=2; i<=n_intrvls_cohort(k); i++)
    {
      //iy   = cal_yr_cohort(k,i); // assign temporary if using multiple times for overhead reasons
      int iy = cal_yr_cohort(k,i)- styr + 1;
      i_intrvl = cal_intrvl_cohort(k,i);
      iage = ages_cohort(k,i);
      F_fishery_cohort(k,i) = fishery_select(iage) * F(iy) * fishery_sample_prob_vec(i_intrvl);
      F_survey_cohort(k,i) = survey_select(iage) * q * survey_sample_prob_vec(i_intrvl) * survey_effort_vec(i_intrvl);
      Z_cohort(k,i) = F_fishery_cohort(k,i) + F_survey_cohort(k,i) + M;
      survival_cohort(k,i) = mfexp(-Z_cohort(k,i)/double(n_intrvls_per_yr));
      cumulative_survival_cohort(k,i) = cumulative_survival_cohort(k,i-1) * survival_cohort(k,i);
      cumulative_nat_mort_cohort(k) += cumulative_survival_cohort(k,i-1)*(1. - survival_cohort(k,i))*M/Z_cohort(k,i);
    }
 }
  ll.initialize();
  for(i = 1; i<=fishery_n_obs; i++)
  {
	  log_exploit = 0.0;
	  ll_i = 0.0;
		  //cout << "ll_i: " << ll_i << endl;
	  if(fishdat_cohort_intrvl_index(i) > 1)
	  	ll_i += log(cumulative_survival_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)-1));
	  log_exploit =  log(1-survival_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));
	  log_exploit += log(F_fishery_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));
	  log_exploit -= log(Z_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));
	  ll_i += log_exploit;
	  ll_i *= fishdat_catch(i);
      ll_i -= gammln(fishdat_catch(i)+1.);
	  ll += ll_i;
  }
  for(i = 1; i<=survey_n_obs; i++)
  {
	  log_exploit = 0.0;
	  ll_i = 0.0;
	  if(survey_cohort_intrvl_index(i) > 1)
	  	ll_i += log(cumulative_survival_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)-1));
	  log_exploit =  log(1-survival_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  log_exploit += log(F_survey_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  log_exploit -= log(Z_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  ll_i += log_exploit;
	  ll_i *= survey_catch(i);
      ll_i -= gammln(survey_catch(i)+1.);
	  ll += ll_i;
  }
  N0 = mfexp(log_N0);
  for (k=1; k<=n_cohorts; k++)
  {
  	ll += gammln(total_fishery_catch_cohort(k) + N0(k) + total_survey_catch_cohort(k) + 1.) - gammln(N0(k)+1.);
  	ll += N0(k)*log(cumulative_survival_cohort(k,n_intrvls_cohort(k)) + cumulative_nat_mort_cohort(k));
  }	
  cout << "ll: " << ll << endl;
  nll = -1.*ll;
}

void model_parameters::report()
{
 adstring ad_tmp=initial_params::get_reportfile_name();
  ofstream report((char*)(adprogram_name + ad_tmp));
  if (!report)
  {
    cerr << "error trying to open report file"  << adprogram_name << ".rep";
    return;
  }
  report << "q" << endl;
  report << q << endl;
  report << "survey_sel50" << endl;
  report << survey_sel50 << endl;
  report << "survey_selslope" << endl;
  report << survey_selslope<< endl;
  report << "fishery_sel50" << endl;
  report << fishery_sel50 << endl;
  report << "fishery_selslope" << endl;
  report << fishery_selslope<< endl;
  report << "M" << endl;
  report << M << endl;
  report <<"N0"<<endl;
  report << N0 + total_fishery_catch_cohort + total_survey_catch_cohort <<endl;
  report <<"F"<<endl;
  report << F <<endl;
}

void model_parameters::preliminary_calculations(void){
  admaster_slave_variable_interface(*this);
}

model_data::~model_data()
{}

model_parameters::~model_parameters()
{}

void model_parameters::final_calcs(void){}

void model_parameters::set_runtime(void){}

#ifdef _BORLANDC_
  extern unsigned _stklen=10000U;
#endif


#ifdef __ZTC__
  extern unsigned int _stack=10000U;
#endif

  long int arrmblsize=0;

int main(int argc,char * argv[])
{
    ad_set_new_handler();
  ad_exit=&ad_boundf;
  gradient_structure::set_MAX_NVAR_OFFSET(1600);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(200000);
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(1800); 
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
  arrmblsize=100000000;
    gradient_structure::set_NO_DERIVATIVES();
    gradient_structure::set_YES_SAVE_VARIABLES_VALUES();
  #if defined(__GNUDOS__) || defined(DOS386) || defined(__DPMI32__)  || \
     defined(__MSVC32__)
      if (!arrmblsize) arrmblsize=150000;
  #else
      if (!arrmblsize) arrmblsize=25000;
  #endif
    model_parameters mp(arrmblsize,argc,argv);
    mp.iprint=10;
    mp.preliminary_calculations();
    mp.computations(argc,argv);
    return 0;
}

extern "C"  {
  void ad_boundf(int i)
  {
    /* so we can stop here */
    exit(i);
  }
}
